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Abstract. We study the thermodynamic properties of a microscopic model of 
coupled oscillators that exhibits a dynamical phase transition from a desynchronized to 
a synchronized phase. We consider two different configurations for the thermodynamic 
forces applied on the oscillators, one resembling the macroscopic power grids, and one 
resembling autonomous molecular motors. We characterize the input and the output 
power as well as the efficiency at maximum power, providing analytic expressions for 
such quantities near the critical coupling strength. We discuss the role of the quenched 
disorder in the thermodynamic force distributions and show that such a disorder may 
lead to an enhancement of the efficiency at maximum power. 
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1. Introduction 

The stochastic thermodynamics of microscopic systems has been the subject of intense 
investigation in recent years [1] in an attempt to extend basic concepts of macroscopic 
classical thermodynamics to the microscopic realm in general, and to out-of-equilibrium 
microscopic systems in particular. Notably a lot of effort has been devoted to the 
characterization of the efficiency of microscopic devices, that can transform heat or 
chemical energy into mechanical work. While for the hrst type of devices (autonomous 
heat engines) the efficiency is bounded by the Carnot limit [2l|3], in the case of isothermal 
engines the efficiency is constrained by the thermodynamic limit 1. In both cases the 
upper limit is reached for quasi static operation, resulting in a vanishing power output. 
Thus, a more relevant quantity to study is the efficiency at maximum power (EMP), that 
exhibits an interesting universal behaviour for different types of devices [SlEllllElinilTllH]. 
In particular the EMP in the linear regime is 1/2 of the maximal allowed value, while 
the behaviour beyond the linear regime depends on the details of the coupling between 
the energy producing and the energy consuming cycles. 

Many of the theoretical studies have been directed toward the characterization of 
the EMP in single devices such as soft nanomachines [5], single molecular motors m, 
devices involving single electron transport [H [9l HO] , or single entropy-driven motors m- 
However, an important class of microscopic devices is represented by cell molecular 
motors, which operate in crowded environments where their mutual interaction can 
become signihcant. For example, many kinesin motors walk on the same microtubule 
leading to traffic jam formation in some case [12], while several motors can pull the same 
cargo resulting in a strong cooperative effect [m m ng usi El]. Furthermore, recent 
studies on synthetic nanomotors have been conducted with the aim of reproducing the 
performance of their biological counterparts [181 EHl 120] . In this regard, it is important 
to note that it is now possible to engineer molecular spiders that exhibit directional 
movement and behave like robots by carrying out a sequence of predetermined actions 
[20] . These artihcial motors might in the future be organized in teams, to optimize, 
for example, their transport properties and efficiency [21]. Thus it is clear that future 
research on thermodynamic property optimization will deal with teams of interacting 
motors, where the dynamical phase these motors operate in becomes relevant. 

One of the most studied models of interacting particles in out-of-equilibrium physics 
is the exclusion process, which exhibits three distinct dynamical phases with different 
densities and particle currents [23]. Furthermore, the exclusion process is often used 
to model molecular motors moving on a lattice, see, e.g. [25l [26]. We have previously 
investigated the issue of EMP in isothermal interacting motors, modelled as an exclusion 
process on a single lattice [8] [22], or on a network [23]. In these studies, we found an 
increase of the EMP in a many-motor system with respect to the single motor case, for 
a suitable choice of the model parameters. Remarkably, in [8] [22] we found that the 
enhancement of the EMP occurs in a range of parameter values compatible with the 
biological estimates for the molecular motor Kinesin. From those studies we concluded 
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that after a dynamical phase transition the dynamical response of the system to an 
external drive can change, leading in turn to a change in the thermodynamic properties. 
Specihcally the dependence of the delivered power on the driving thermodynamic forces 
may vary. 

One of the main limitation that one faces when studying the thermodynamic 
properties of exclusion processes on a lattice, is that the intensity of interaction between 
the motors can only be indirectly tuned by changing the kinetic parameters and thus the 
density of motors on the lattice P1221123]. Furthermore, if one wants to study the effect 
of force disorder on the motor particles, one has to resort to numerical simulations, as 
no exact result exists for the exclusion process with heterogeneous particles. 

Instead, here we consider a model of N interacting microscopic particles, where the 
particle-particle interaction is an explicit parameter, that can be tuned in order to drive a 
dynamical phase transition, from a weakly interacting - incoherent system to a strongly 
interacting - coherent system. This model was originally introduced by Sakaguchi in 
izg as an extension of the Kuramoto model (KM) j2H], to study the synchronization of 
a group of interacting oscillators in contact with a reservoir at constant temperature. 
Furthermore, the effect of quenched disorder in the thermodynamic force distribution 
can be taken into account within the present model. The model is introduced and 
discussed in section |2l The N interacting particles can be viewed as a network of energy 
producers and users or as a system of interacting autonomous motors under the effect 
of thermodynamic forces. Since the dynamical phase diagram can be obtained in terms 
of the particle interaction strength, temperature and force distribution, in section |4] 
we will discuss how to calculate the relevant thermodynamic quantities, namely the 
deliverer and input power, and the efficiency. We will consider two possible scenarios 
as far as the force distribution is concerned. In the hrst one, sec. 14.11 either a positive 
or a negative force is applied on each particle. In the second scenario, sec. 14.21 both a 
positive and a negative force is applied on the same particle. We will thus discuss how to 
optimize the delivered power for the different types of network models, hence obtaining 
the EMP in terms of the interaction intensity, and thus of the coherence between the 
particles’ motion. We will hnally discuss the effect of the quenched disorder in the force 
distribution on the thermodynamic quantities. 

Interestingly, the model that we use here is a microscopic version of a model used to 
mimic macroscopic power grids. Indeed the dynamics of interconnected power grids can 
be mathematically represented by a complex network of coupled oscillators [291130] . while 
at the macroscopic level one faces optimization problems different from the microscopic 
case, as shortly discussed in section |3l 
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2. The Sakaguchi model 


We consider a system of N coupled oscillators, originally introduced by Sakaguchi en, 
described by the Langevin equation 

= fi- - 0i(i)) + Viit), (1) 

j 

where ft is an external constant force, and the Gaussian noise rji obeys the fluctuation- 
dissipation relation 

{Vi{t)Vj{t')) = 2kBT6ij6{t - t'). (2) 

Notice that we have chosen the system units such that the external force fi has dimension 
of frequency, which corresponds to taking the friction coefficient in eq. ([2]) equal to one. 
By introducing the complex order parameter 

a(t)exp(i^/^(t)) = ^^exp(i0j(t)), (3) 

j 

where 0 < a{t) < 1 measures the system coherence and is the common average 
phase, equation ([T]) becomes 

^i = fi- Ka sin(0i - V’) + Vh (4) 

where we understood the dependence on time. Let /o be the mean deterministic force, 
calculated over the N oscillator sample /o = J2jfj/^y expect that the center of 
mass will oscillate with the frequency /o, so we can set 'if(t) = fot + -00 and thus we can 
redehne the dynamical variables as 9i = (fi — 'ifif), so as eq. (jlj) reads 

6i = Ui- Ka sin(6'j) + rji, (5) 

where we have redehned the external force as (^i = fi — fo- 

In principle eq. ([5]) represents a set of N coupled equations for the variables 9i, 
since a is given by eq. ([3]). However, as iV —)■ oo, one can replace the actual value of a 
with its mean held value, and so eq. ([5]) becomes uncoupled. Such a mean held value 
can be obtained self consistently as discussed below. Eq. ([5]) corresponds to a Brownian 
particle moving in a periodic potential under the ehect of a constant drift force cjj. Here 
and in the following we assume that the system reaches a steady state in the long time 
limit. In the course of this paper, we will discuss this assumption where relevant. The 
Langevin equation can be reformulated in terms of a Fokker-Planck (FP) equation for 
the probability distribution function (PDF) of hnding the particle i at position 9 at time 
t 


dtp{9,uji,t) = de [(iPcrsin^ - oJi)p + Td0p]. 


( 6 ) 


Thus, the stationary probability distribution function (PDF) of the position of such a 
particle reads nm 


p{9,uji) 


/(27r) 

1 — exp {—(327iUi) 





( 7 ) 
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where I{x) = dyexp[(3{Kacosy + Uii/)], and A/" is a normalization constant 
depending implicitly on (3 = 1/T, K ■ a and cuj. The steady-state probability current 
thus reads Jss = A/", and the particle steady-state velocity reads 

V0{Ka,uJi,T) = 2 ttM 


f27r 


= 2t;{I3 / 


I{2n) 


1 — exp {—(32TTUJi) 


-m 


( 8 ) 


As A^ —)■ oo, we can adopt a continuous description, where the constant forces 
acting on the oscillators are distributed according to the probability distribution g{f) 
with mean value /q. By introducing the shifted force distribution 

= 5'(/o + 1^); (9) 

the self-consistent equation for the modulus a of the complex order parameter, 
characterizing the degree of order or coherence in the conhguration of the variables 
6i is then given by 


r-27r 


ere 


= dugoiu) / d0p(0,a;)exp(i0) 


which can be decomposed into its real and imaginary part 


r. 27 r 


cr cos('0o) = / 9 o{^) / d6*p(6*,a;) cos^ , 


f27r 


cr sin('0o) = / du go{oj) / d6*p(6*,a;) sin6*. 


( 10 ) 


( 11 ) 


( 12 ) 


By assuming that the force distribution g{f) is symmetric around /o, and noticing 
that p{6, —u) = p{—0, u), the imaginary part on the right-hand side of eq. ffTOj) vanishes, 
and so one is left with 


r2n 


cr = / da;5fo(i^) / d9p{9,u) cos{9), 


(13) 


whose solution provides the mean held value for a. 

As discussed in [2^, for N ^ oo this model exhibits a critical coupling strength 


Kc, such that for K > Kc the systems exhibits a dynamical phase transition with 
synchronization a > 0, while the system is incoherent for K < Kc-, and each particle 
described by the coordinate 9i oscillates with its proper frequency oji. Thus, for K > Kc 
we expect a to be positive but small, and we can expand eq. flT^ in powers of e = Ka/T, 
obtaining 


KaT 


b + CX) 


a = 




L (T^ + 


1 - 


(T2 - 2u^) 


2 (T2 + a;2) (dT^ + a;^) 

(32-4 _ 172-2^2 ^ 4^4) n 

+ 0 (e") , 


while expanding eq. (jH]) 

Vg{u) = U 


"^4 (T2 + (4T2 + u;2) (QT^ + . 

the average velocity of the dynamical variable 9 reads 

AV^ (5T2 - luY 
~ 2(T2 + a;2) ^ 8 {T^ + (4T2 + 

(23T^ - 24T2a;2 + o;^) 


(14) 


16 (T2 + uY (4T2 + uY (9T2 + . 


+ 0(e«). 


(15) 
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It is worth to note that the hrst two coefficients of the expansion fll4l) were recently re¬ 
derived in ref. [33] by using a different approach: the author mapped the deterministic 
time evolution of the Kuramoto model order parameter into a stochastic process as 
given by eq. o. and applied a fluctuation relation to the thermodynamic irreversible 
work done on such a stochastic system. 

Inspection of eq. (na provides the critical coupling strength for which a non¬ 
vanishing solution to that equation appears 



(16) 


The value of the order parameter a as a function of K and T, for K > can be 
obtained by solving eq. (fT4|) . which gives 



where we have introduced the quantity AK = K — Kc, and 



(18) 


(19) 


The hgures contained in this manuscript and discussed in the following, are made 
by using the approximate expressions (na and 031. and choosing the values of the 
parameters such that a < 0.5. 

A few comments on the velocity vg{uj) are now in order. Such a quantity represents 
the velocity of the dynamical variable 6 and is thus the velocity deviation of a particle 
under the effect of a force uj with respect to the center of mass velocity /q. Inspection 
of eq. (I5|) suggests that for K < Kc (thus for a = 0) vg{u) = oo. On the other hand 
Vg{u}) goes to zero as K increases above Kc- the higher K the higher are the barriers of 
the periodic force in eq. (]5|), while a is also an increasing function of K. 

3. Macroscopic power grids 

Here we want to establish contact between the network of microscopic oscillators 
discussed in the previous section, and a corresponding macroscopic model used to 
describe macroscopic power grids. This analogy will be useful in the next section, where 
we will introduce the thermodynamic forces and powers characterizing the microscopic 
model. Given that in the macroscopic realm, one deals with alternating current (AC) 
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networks, in terms of equations, this amounts to write a system of equations for the phase 
angle (pi for both the generators and the users [30] which correspond to an extended KM 

N 

Mipi + Dipi = uJi — ttij sin((/)i — <pj); if Ms a generator (20) 

i=i 

N 

Dipi = oji — ttij sin(0j — pj), if Ms a consumer; (21) 

with inertia coefficients Mj (representing, e.g., the large rotational inertia in turbine 
generators), viscous damping Di, power injection Ui (consumption if Ui < 0) and power 
flows along the lines aijsm{pi — pj) with coupling strength aij. This model exhibits a 
synchronized phase that corresponds to a power grid that operates in a steady state 
with spatially uniform frequency. 

In this kind of power network one may want to determine the optimal operation 
conditions to, e.g., avoid blackouts. For example in a static grid with a fixed number 
of producers and users, where the generators have a maximum power capability, and 
the users are characterized by a well known average consumption, one may want to 
optimize the grid in order to avoid that the consumers’ load exceed the generators 
capability, leading to blackouts. This can be done by using optimization algorithms for 
graphs, see, e.g., [34] . 

Another, perhaps more interesting optimization problem considers a power grid 
as a dynamical system, where both energy producers and users can be dynamically 
connected or disconnected over a given time period. On the producers side, this is the 
case of renewable energy sources, such as wind and solar power, which are stochastic 
in nature and often uncontrollable, resulting in severe difficulties in the maintenance 
of the balance between load and generation [35]. Thus in the future opportunistic 
users can access the energy system according to the availability of system resources and 
differently from the ” always-on” demand of traditional energy users, their consumption 
can exhibit peaks of activity. In this scenario, the challenge is to coordinate and manage 
dynamically interacting power grid participants. 

4. Stochastic Thermodynamics of the microscopic model 

In section [2] we have discussed the dynamical properties of the model system we will use 
in the present paper. We can now turn our attention to its thermodynamic properties, 
namely the input and delivered power, and the system efficiency as a global motor. 

One can consider two possible scenarios as far as the forces applied on each single 
particle are concerned. 

In the hrst case the forces acting on the motor system can be either positive or 
negative, thus resembling the macroscopic power grids of power plants and consumers 
considered, e.g., in [361 133 ED] EE]- Differently from those works, we consider here 
microscopic oscillators, in the over-damped regime, and with white noise acting on 
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them. In this scenario, taking inspiration from the macroscopic realm, one may call 


users those oscillators with a negative force acting on them /j < 0, and producers those 


oscillators with a positive force /* > 0, and a single force distribution characterizes the 
system. 

The second possible scenario resembles the case of molecular motors, where both 
a negative < 0) and a positive force > 0) are applied on the same particle i. 
This is the case in, e.g., biological molecular motors such as kinesin and myosin |39l 00] 
where the energy extracted by ATP hydrolysis drives the motor forward (corresponding 
to > 0) while the motor does work to carry a cargo, modelled by a negative load 
(corresponding to f~ < 0) In this case one deals with two different distributions of 
forces, g+{f+) and (;_(/_). The homogeneous case, where the same positive /+ and 
negative force /_ where applied on all the motors, modelled as diffusing particle on a 
lattice with an exclusion rule, was studied, e.g., in [8] l2^ 123]. 

In both scenarios, in order for the system to behave globally as a motor, and to 
perform work against the negative forces, we must require the center of mass to have 
an average positive velocity, and thus /o > 0. 

In the following, we will consider, for both scenarios, the delivered power Pout 
and the input power Pin, and optimize Pout wrt different parameters. We will 
thus characterize the efficiency at maximum power (EMP) rj* = Pout/^in Plainly in 
proximity of the dynamical phase transition, and where possible, we will discuss the 
thermodynamic properties of the system in the whole range of parameter space. 

4.I. Single force distribution 

We consider here the stochastic thermodynamics of a system with either negative or 
positive forces applied on each oscillators, and distributed according to the single PDF 


aif)- 


We can thus introduce the relevant thermodynamic quantities, namely the average 
input power, absorbed by the producers, and the average output power released by 
the users. Recalling that ve as given by eq. ([HD gives the deviation of the i-th particle’s 
average velocity from the center of mass velocity /o, the average output and input power 
read 



( 22 ) 


( 23 ) 
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while the thermodynamic efficiency of the system reads 


V 


P, 


out 


Pw 


(24) 


Substituting eqs. (IT^ and flTTD into (12^ and fl25]) the output and input power 
becomes, up to the second order in 


where 


„ , K-K,,, , (K - IQ-'- (Ifh + hQ - QI'ilQ 

— -Pf\ ~r -TT^TT--i o “T - 




mt 


r, r,> , - !<■= ,> , (A' - KQ- (I>h + hQ - Ql'ilQ) 

■'in — -^0 X -'9 “T 




IlKt 


p< — 
Aq — 


12 = 


I4 = 


r-fo 

- dugo{u){cj + fof <0, 

J —OO 

f A z' ^ ^ n 

/ <^^9o{^)iT77pr—p^ > 0 , 

' —OO 

r-fo 


dujgoiu 


2(T2 + a;2) 

{u + fo)u (ST^ - o;^) 


8{T^+u^Y{4:T^ + uj^y 


(25) 

(26) 

(27) 

(28) 
(29) 


with analogous definitions for P^, and . We notice that, in absence of partial 
synchronization {K < Kc, a = 0), i.e., when the users and the producers are decoupled, 
eq. fl22l) . fl25|) . and fimi predict that the delivered power is negative. Since vo{u) = to 
for K < Kc, the users oscillates with their proper frequency (force) which is negative, 
and so the product of the applied forces times the average velocity is positive. The term 
vq{uj) in eq. (122|) is always negative, as the integration variable runs over negative value. 
Thus by increasing K above Kc the modulus of V 0 {(jj) decreases and tends to zero for 
very large K. This implies that for some value of K the rhs of eq. fl2^ becomes positive, 
such values depending on /o and T, and on the details of the distribution goiu), e.g. its 
width. 


4-1.1. Optimization Here we aim at optimizing the delivered power eq. fl2^ wrt some 
of the relevant parameters. 

Optimization wrt the coupling strength dPont/dK = 0 gives: 


dp 


out 


"-/o 


dK 


dcj go{u)udKVe{K, uj,T). 


(30) 


Recalling that the average velocity ve goes to zero as K increases above K^ we find 
that, for (u < 0, Ve{u},K) is an increasing function of K, ranging from oo for K < Kc 
and approaching zero as K ^ oo. Thus, from eq. (1301) it follows 


dPont _ f if R' < Kc, 

dK ~ \ > 0 if R: > Kc. 
Similarly one finds 


(31) 


dP{a _ R’ liK < Kc. 
Hk ~ I < 0 if A' > Kc . 


( 32 ) 
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Thus, if K is the free parameter, the optimal delivered power is achieved for 
K ^ oo, corresponding to the limit of strong coupling between users and producers, 
with an EMP rj* = (/_) / (/+) as obtained by eqs. (j2^ and (j23|) . where (/_) and (/+) 
are the average negative and positive forces, respectively. 

No similar inequalities can be found when one tries to maximize Pout with respect 
to other parameters, for example /q. So one should consider specihc cases for the force 
distribution in order to study the relevant thermodynamic quantities. 


4-1.2. A specific distribution In order to exemplify the results discussed in this section, 
here we consider the specihc distribution 

9{f) = ^ [S{f - (/o + s)) + d{f - (/o - s))], (33) 

where is the variance of the distribution, with s > /o > 0, i.e., there are just two 
types of oscillator, the users with an applied force /o — s < 0 and the producers with an 
applied force /o + s > 0. The shifted force distribution thus reads 

fi'o(^) = ^ [<5(cu - s) + 5(a; + s)]. (34) 

For such a distribution the critical coupling strength reads 


Pc = 2 


(S2 + T2) 

T 


(35) 


This corresponds to the bimodal distribution considered in HU, where the linear 
stability of the incoherent solution p{0,u)) = l/27r of the FP equation ([6]) was studied, 
corresponding to the non-synchronized phase a = 0. The authors showed that in the 
limit N ^ oo the FP equation eq. ([6]) exhibits a steady state solution for s < T and 
K > Kc, while for s > T and K > AT the system exhibits an oscillatory state. In 
the following we will thus take s < T, and use the steady state solution ([7]) to eq. (E]). 
Equations fl25|) - fl26|) thus become 


Pont = ^is - fo){fo-ve{s)), (36) 

Pin = +/o)(/o +'i^e(s)), (37) 


We recall that for K < Kc (i.e. for a = 0), vq{s) = s, and because of the condition 
s > /o we have Pout < 0, i.e. when the producers and users are not coupled, the users 
oscillates with their proper frequency fo — s resulting in a negative Pout • The delivered 
power will become positive for some value of P > Pc, when vo{s) in eq. fl36|) becomes 
smaller than /q. 


4.1.3. Optimization Here we optimize the delivered power for the force distribution 
(l3T)) . which corresponds to a system where on each oscillator there is either a positive 
fo + soia, negative fo — s force with probability 1/2. From eq. flTbl) we easily obtain 
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the expression for the velocity deviation from the center of mass up to the fourth order 
in a 


(ST^ - s^) 

2(s2 + T2) ^ 8 (s2 + T2)^(s2 + 4T2) 

A'®a® ( 23 T^ - 24 T 2 s 2 + 54^ 

16 (T2 + s 2)3 (4T2 ^ (-9^2 + 52) J ’ 


while the order parameter, as given by eq. (Ell, becomes 


a = 


'AA'T(s2 + 4T2) 
K^{T^ - 2s2) 


1 + AK 


8s® + 97- dOs^T^ + 15T®' 
2Kc{T^ -2s^y{s^ + 9T^) \ ’ 


(39) 


up to the order 3/2 in AK. 

We can now optimize Pout, as given by eq. fl36|) . wrt to different parameters: 
i) By optimizing wrt to K at hxed /o and s: OkPoui = 0, one obtains 


dPr 


out 


dK 


- 1 / 2 ( 5 -/o) 


dv0{s) 

dK 


> 0 , 


(40) 


since ve^s) is a decreasing function of K, as already discussed above for a general force 
distribution go{u)). 

ii) By optimizing Pout wrt the average force, at hxed K and s, df^Pont = 0, one 
obtains 


as.K) 


ve{,s) + s 
2 


(41) 


and since the condition s > vg{s) > 0 holds for any K > 0, we have s > /q (s, K) > s/2. 



Figure 1. EMP ij* as obtained by maximizing Pout wrt /o, as a function of K for 
different values of the quenched disorder standard deviation s. Here T = 1. The 
dashed lines correspond to the approximated expression (1451) . 


We can thus calculate the delivered and the input power, and the efhciency at the 
maximum 

P„, = -«»(«))" 


(42) 
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K = + 'y0(s))(s + 3u0(s)) 

, ^ _ AK^ {s^ + 4T")^ 

^ (3s + t>e(s))(s + 3t>0(s)) 16K^{T^-2s‘^f 

^AK‘^1'1 5 s 2 \ 

- r2 Vi ^ ’ 


(43) 

(44) 

(45) 


where we have used (|3^ and (13^ to expand t]* up to the lowest order in AK and s/T. 
Plots of rj* as a function of K for different values of s are shown in £g. [1] Inspection 
of this figure, as well as of eqs. fl42l) . fl43ll and fj44ll suggests that, for fixed s, when K 
increases above K^, the optimal output power fl42|) increases, the optimal input power 
fl43|) decreases, and this results in an increase of the EMP (jH]). This is a consequence 
of the fact that ve{s) —?■ 0 in the limit K —)■ cxo, where rj* = 1/3. 

Inspection of hgure [H as well as of eq. fl45|) . suggests that a higher degree of 
quenched disorder, as parametrized by s, leads to a larger EMP close to the dynamical 
phase transition. However, the analysis of the behaviour of eqs. fl42|) . fl43|) and fl44)) 
at fixed AK and varying s is not so straightforward. Graphical analysis of eqs. (I42ll . 
(|43|) (not shown) indicates that both and P,* increase with s, with increasing 
faster. This graphical check can be done in the range of parameters where eqs. (1381) 
and fl3^ holds, i.e. close to the critical point. However, it is worth to note that for a 
hxed AK, one hnds P/ut('S = 0) = 0. Furthermore inspection of eq. (j5]) also suggests 
that vg{s) —)■ s, as s —)■ oo, and being P/(,^ a positive quantity, it must have at least one 
maximum for s G [0, +oo[. On the other hand, P,* , as given by eq. (I43|) in an increasing 
function of s. Accordingly rj* has at least one maximum for s G [0, +cxd[. 

Hi) In order to optimize Pout, as given by eq. fl36l) . with respect to the quenched 
disorder standard deviation s, one has to solve the equation dgPout = 0. This equation 
has no analytic solution s*, but it can be solved numerically, in order to find the EMP 
for different values of /o and K, as shown in fig. |2l Still one has the constraint s* > /o, 
in order for the force on the user to be negative. 

One can also hnd an approximate expression for s* at the lowest order in AK and 
s/T as follows. 

By solving dgPout = 0 for K < Kc, one easily finds that the maximum is given by 
■Sq = /o, cind thus P/^^ = 0. For K > K^, one can derive -Pout wrt s, and then expand 
the equations up to the second order in {K — Kc) and As = s — Sg, solving for s, and 
plugging the value s* that maximize Pout into the expression for r], one finds 


{K - Kc)^T^ (4T^ + f^)^ (T^ + 7f^) 

64 (T2 - 2/2) (T2 + /2) (T4 + 3T2/o' - f^)" 
(I 5/2 \ 

T2 V4 8T2y ’ 


(46) 


where the last expression gives r]* to the lowest order in AK and fo/T. It is worth 
noting that the coefficients in the series expansions are the same as in eq. fl45]) . 

Note that AK depends implicitly on s* through Kc, as given by eq. (13^ . when one 
replaces s with s*. Inspection of eq. (15B]) suggests that the optimal quenched disorder 
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Figure 2. EMP rj* as obtained by maximizing Pout wrt s, as a function of /o for 
different values of the coupling constant K. The dashed lines correspond to the 
approximated expression (1461) . 


standard deviation s* increases as /o increases, but this in turn leads to a smaller value 
of AK for fixed K, as Kc{s = s*) also increases. Thus increasing /o, and optimizing 
-Pout wrt s drives the system towards smaller values of synchronization a, resulting in 
a smaller rj*. This is in agreement with the results reported in £g. [2l showing that the 
EMP close to the dynamical phase transition, is enhanced by decreasing the applied 
average force /q. 


4-1.4- Gaussian distribution We now consider the following distribution for the force 

(47) 


9(1^) = 


e 2 s^ 


\/2 


TTS^ 


From eq. (1T6|) we can calculate the critical coupling strength 


Kr = 2\ —s 


e ^ 


“f I A. 


(48) 


and from eqs. and we can calculate the output and the input power for K < Kc 
which read 


with 


Pout,0 ®/o 


Pin,0 — s/o 


_ fo 
e ^ 

f2 

e ^ 


2 


+ fo) 
+ fo) 


1 — erf 


1 + erf 


fo 

y/2s 

fo 

y/2s 


hniPout,o = 0, 


s —>0 


limPin,o =/o, 




2 ’ 


<0, 

(49) 

>0, 

(50) 


(51) 


(52) 


2 
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No analytic result can be obtained for Pout, and Pin and thus for the EMP when K > Kc 
for the Gaussian force distribution fl47|) . at variance with what has been done in 
the previous section. However, one can resort to numerical calculations, to integrate 
numerically eqs. (122|) - (|23|) . find the optimal value of Pout wrt to some parameter, and 
thus calculate the EMP. This procedure has been followed in order to calculate the EMP 
as obtained by maximizing Pout wrt to the mean force /q. The results are reported in 
fig. O while for small AiP the EMP is larger the smaller the variance, for sufficiently 
large AJP we hnd the same tendency observed for the bimodal distribution, namely rj* 
as a function of AK is larger, the broader is the force distribution, and thus the degree 
of quenched disorder. 



Figure 3. EMP ij*, as obtained by maximizing Pout wrt /o, as a function of 
{K — Kc)/Kc, for the Gaussian force distribution (HTl) for different values of the 
variance. 


4-2. Distribution of positive and negative forces on the same particle 

In this section, we consider the case where on the same particle, whose dynamics is 
described by eq. two forces are applied, /* _ < 0 and > 0, distributed with two 
PDFs g-{f-) and g+{f+). In this framework the output and input power read 

^+oo fO 

Pout = - / d/+g+(/+) / d/_g_(/_) /_ IvsiU + /- - /o) + A] .(53) 

Jo J —oo 

P + CXD pO 

Pin = / d/+9+(/+) / d/_9_(/_)/+|t,»(/++ /_-/„) + /„]. (54) 

^0 J —oo 

respectively. 

4.2.1. No disorder We start our analysis by considering the trivial case where the same 
forces /+ and /_ are applied on all the particles, with /o = /+ + /_> 0 and /_ < 0. 
So, the total force distribution reads giu) = 6{u) — fo), and the delivered power (l53|) 
becomes 

Pout = ~ 


/_./„ = -/_(/_ + /+), 


(55) 
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which is independent of the coupling K. Thus by maximizing -Pout wrt /_ one finds that 
the EMP is always rj* = 1/2, as in the linear regime case [21 |6l 0 02]. 


4-2.2. Bimodal negative force distribution In order to increase the complexity of our 
system, we consider here the following distributions 


9-^(/) = MJ - /+), 


so as the total force distribution reads 

9(/) = 11'5(/ - (/+ + h-))+Hf- (/+ -t h-)] ■ 


We introduce the following variables 


/l - + h- ^ n 
X = —^^ < 0, 


y = 


h- - / 2 ,- 


2 ’ " 2 
Given the expression for the average force, we obtain a condition on x: 

fo = f++ t = f+ + x >0 => -/+ < a; < 0. 


The delivered power eq. fl53|) thus becomes 
Pont = - [x{x + /+) + yve (y)], 
while the input power reads 

Pl„ = f+{x + f+). 


(56) 

(57) 

(58) 

( 69 ) 

(60) 

(61) 

(62) 

(63) 


It turns out that with this choice of the force distributions the order parameter a and 
thus V 0 depends only on the negative forces, indeed we have 


9 o{(^) = 9 {^ + /o) = 2 ['^(^ “ (/i.- “ /2,-)/2) + - (/2,- - /i,-)/2)] 


= 2 -y) + + y)] 1 


(64) 


and goioj) is the function appearing in the self consistent eq. (ITU]) . 

By deriving Pout with respect to the disorder degree parameter y, we obtain 
dyPont = —{ve {y)+yPe (y)) < where we have assumed y > 0 without loss of generality, 
i.e. Pout decreases monotonically with y. Thus the optimal Pout is trivially obtained for 
y = 0 which corresponds to the case with no disorder, with rj* = /_//+. 

We now optimize Pout eq. dUU]) wrt the variable x, which is equal to the mean 
negative force, and find 

* - _l± 

^ 2 ’ 

and thus 

1 


(65) 


P/ut = 4/+ “ (y ), 


f2 

p* = 

in 2 ’ 


( 66 ) 
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and finally the EMP reads 

fl - ^yve {y) 


V = 




(67) 


We notice that while for K < one finds veiy) = ?/, and therefore in the nnconpled 
regime rj* < 1/2, for K > Kc the velocity V 0 {y) is a decreasing function of K, thus for 
a fixed y the coupling reduces the spread around the mean velocity /o, and thus in the 
limit of large K one recovers the single particle EMP rj* = 1/2. Similarly, (|66l) is an 
increasing function of K, because of the decreasing behaviour of V 0 {y). By expanding 
vg in powers of iP, and noticing that the order parameter a is given by eq. fl3^ , with 
the substitution s ^ y, we obtain 


T] ~ 


fl - 4:y^ , 4 |/ 


+ 


rAJP, 


2fl ■ Tfi-- 

to the lowest order in AJP and y. We notice that the quantity y"^ is equal to the 
variance of the negative forces’ distribution, thus the introduction of disorder in the 
force distribution, on the one hand reduces the EMP for an uncoupled system [K < 
hrst term on the rhs of eq. fl65]) i. on the other hand it increases the slope of the EMP 
above the critical coupling. 

In Eg. mthe EMP rj* is plotted as a function of K for different values of y. 



Figure 4. EMP as a function of the coupling strength K as obtained by maximizing 
Pont wrt the average negative force x, with T = 1, and /+ = 1. The dashed lines 
correspond to the approximate expression (1681) . 


4-2.3. General case: optimization We now consider the case where the system exhibits 
a distribution of both positive and negative forces g-{f-) gVf+)- The distribution 
of the total forces on each particle thus reads 

g(/) = j d/_ d/+ 9 _(/_)g+(/jy/ - (/+ + /_)), (69) 

with 


/o 



d/_ (/_(/_)/_ + 


dUgVU) = f- + U, 


(70) 
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and 


soM = 9(11'+ A) = / <if-<if^-s-(f-)aAf+)s(<^-(y+ + y-)), 


(71) 


where y± = f± — f±. The distribution go{uj) is symmetric around cu = 0 if both the 
distributions g±{f±) are symmetric around their respective average value f±, a symmetry 
that we assume in the following. 

We have thus 


b + CX) 


-Pout — ~ 


dy+g+{f^ 


= - /_ (/+ + /_) 


J —00 
^+00 

'-/+ 


dy-g-{f-) f- [ve{f+ + /- - /o) + /o] 
r-f- 


dy+g+{y^ 


and similarly 


Pin — /+ {f+ + /-) + 


b+cx) 


'-/+ 


dy+g+{y^ 


a-f- 


dy-g-{f-) y-ve{y+ + y-), 


dy-g-{y-) y+ve{y+ + y-), 


where in the last equality we have used the above mentioned symmetry oig±{f±). Thus, 
if we want to optimize Pout wrt to /_, we obtain 

/+ 


and 


r- = -^- 


f2 

P* _ J + _ 
■^out ^ 

f2 

P* = _|_ 

in 2 ^ 


r+00 

'-/+ 

c 


dy+g+iy^ 


dy+g+{y+) 


r 

' —00 

/■-/- 


dy_g_(y_)y_ve(y+ + y-), 


iy-g-{y-)y+vf{y+ + y-). 


(72) 


(73) 


(74) 


We can thus evaluate the optimal delivered power below the critical coupling, and for 
very large coupling constant 


P* _ 

■^out 


f2 

fl 
4 ’ 


for K < Kc 
for P > Pe 


(75) 


Close to the critical K up to the second order in a, by using eq. 03. we hnd 
Pout = ^ - (y-) + / dy+dy_g+{y+)g_{y_)-'^~^'^~ ^ ^ 


2[T^ + {y_ + y+y] 


fl ^ 


, /2AP 

{yl) -1 


T 


(76) 


p* ^l M / ^\ M f A A ( ^ \ y+(y-+ 

Pin = w + \y+/ + / dy+dy_g+{y+)g_{y_)-^ 


2 

fl 

2 


2[T^ + {y_ + y+y] 


{yl)[ 


f2AK 


T 


+ 1 


(77) 
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to the lowest order in AK/T. We obtain finally the EMP 

. ^ , 2AKfl{2(£)-{y'i)) + S{y'i)(yl) 

' T {fl + 2{yl})^ ■ ' ' 

inspection of this last equation suggests that the maximal slope of rj* as a function of 
AK is obtained for (|/+) = 0. As far as the variance of the negative forces is concerned, 
we find a similar scenario as in the previous section: while (|/i) reduces the EMP for 
the uncoupled system, above the critical coupling the EMP increases faster the larger 

is (y-)- 

5. Conclusions 

In the present paper we have investigated the thermodynamic properties of a model of 
microscopic oscillators, subject to thermodynamic forces. We considered the effect of 
the disorder on the delivered and injected power and on the EMP, and discussed the 
critical behavior of such quantities for different force distributions. 

We considered two forces distribution types, one that resembles the macroscopic 
power grids, and one that resembles a system of interacting autonomous motors. 

For the first type of force distribution we find that, at fixed coupling strength, a 
larger degree of disorder leads to an increase in the EMP, at least close to the critical 
point. 

For the second type of force distribution we find that while the disorder reduces 
both the optimal Pout and the EMP below the critical coupling, above the critical point 
the EMP rate as a function of AK increases as the degree of disorder increases. 

Thus, ideally the system with the optimal thermodynamic performances is 
characterized by a strong coupling {K oo) or absence of force disorder. However, 
in a real system one may have to deal with a finite coupling strength and an intrinsic 
non-vanishing degree of disorder in the force distribution. The results contained in this 
paper characterize the thermodynamic properties of such systems. 

While we were able to calculate the expansion of the energy rates and of the EMP 
close to the critical point, we found that the EMP does not exhibit any universal 
behaviour, at variance with the single device case. On the contrary, the results contained 
in this paper depends strongly on the details of the force distribution. For example, in 
eq. (j78j) . one recovers the EMP value in the linear regime [r]* = 1/2) only when the 
disorder vanishes. 

One of the limitations of the present model is that it exhibits an ”all-to-all” 
coupling, while in a real system the interaction can depend on the distance between 
the network nodes and on the network topology, thus one should replace the interaction 
strength K in eq. ([I]) with an interaction matrix Kij. The thermodynamic properties of 
this extended model are certainly worth to investigate. 

Furthermore, the values of the EMP reported in figs.[Tl[2l[3]are quite small, thus the 
characterization of the response of the network injected and delivered power, or of its 
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efficiency to a change in the network topology is certainly worthy of future investigation. 
For example, one may want to hnd the connectivity matrix between the different nodes 
that optimize the relevant thermodynamic quantities. 
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